Changing the Spatial Extent or resolution of rasters

Cropping Rasters

  • Cropping is the production of a smaller raster from an existing one by selecting a certain range of rows and columns.

  • Subsetting is also a form of cropping,although in practice we select the required range of rows and columns by overlaying the raster with another spatial layer.

  • The crop function, given a raster and an Extent object, returns a smaller raster with the required extent.

  • Instead of an Extent object, a raster or a vector layer can also be provided, in which case their extent is extracted and used in cropping.

We’ll show this by an example but first we’ll import the required libraries. 

library(raster)
library(sp)
library(rgdal)
library(dplyr)
library(magrittr)

Next, we’ll read in a raster of Bangalore.

raster_file_list <- list.files(path = "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1/", pattern = "[0-9].tif$",ignore.case = TRUE, full.names = TRUE) 

raster_file_list
##  [1] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B1.TIF" 
##  [2] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B10.TIF"
##  [3] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B11.TIF"
##  [4] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B2.TIF" 
##  [5] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B3.TIF" 
##  [6] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B4.TIF" 
##  [7] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B5.TIF" 
##  [8] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B6.TIF" 
##  [9] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B7.TIF" 
## [10] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B8.TIF" 
## [11] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B9.TIF"

We don’t want band 8 for this analysis since the resolution of band 8 (15m) is different from other bands (30m)

raster_file_list <- raster_file_list %>% grep (pattern = "B8", x = ., ignore.case = TRUE, invert = TRUE, value = TRUE)
raster_file_list
##  [1] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B1.TIF" 
##  [2] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B10.TIF"
##  [3] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B11.TIF"
##  [4] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B2.TIF" 
##  [5] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B3.TIF" 
##  [6] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B4.TIF" 
##  [7] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B5.TIF" 
##  [8] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B6.TIF" 
##  [9] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B7.TIF" 
## [10] "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1//LC08_L1TP_144051_20181020_20181031_01_T1_B9.TIF"
rasterStack <- stack(raster_file_list)
rasterStack
## class       : RasterStack 
## dimensions  : 7711, 7551, 58225761, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 632985, 859515, 1323585, 1554915  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names       : LC08_L1TP//1_01_T1_B1, LC08_L1TP//_01_T1_B10, LC08_L1TP//_01_T1_B11, LC08_L1TP//1_01_T1_B2, LC08_L1TP//1_01_T1_B3, LC08_L1TP//1_01_T1_B4, LC08_L1TP//1_01_T1_B5, LC08_L1TP//1_01_T1_B6, LC08_L1TP//1_01_T1_B7, LC08_L1TP//1_01_T1_B9 
## min values  :                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0 
## max values  :                 65535,                 65535,                 65535,                 65535,                 65535,                 65535,                 65535,                 65535,                 65535,                 65535
rasterBrick <- brick("/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1/rasterBrick_b1_to_b11.grd")
rasterBrick
## class       : RasterBrick 
## dimensions  : 7711, 7551, 58225761, 10  (nrow, ncol, ncell, nlayers)
## resolution  : 30, 30  (x, y)
## extent      : 632985, 859515, 1323585, 1554915  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1/rasterBrick_b1_to_b11.grd 
## names       : LC08_L1TP//1_01_T1_B1, LC08_L1TP//_01_T1_B10, LC08_L1TP//_01_T1_B11, LC08_L1TP//1_01_T1_B2, LC08_L1TP//1_01_T1_B3, LC08_L1TP//1_01_T1_B4, LC08_L1TP//1_01_T1_B5, LC08_L1TP//1_01_T1_B6, LC08_L1TP//1_01_T1_B7, LC08_L1TP//1_01_T1_B9 
## min values  :                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0,                     0 
## max values  :                 41084,                 31902,                 27503,                 50157,                 51469,                 59522,                 65535,                 65535,                 65535,                 13633
names(rasterBrick) <- gsub(pattern = "LC08_L1TP_144051_20181020_20181031_01_T1_B", replacement = "band_", x = names(rasterBrick), ignore.case = TRUE)
names(rasterBrick)
##  [1] "band_1"  "band_10" "band_11" "band_2"  "band_3"  "band_4"  "band_5" 
##  [8] "band_6"  "band_7"  "band_9"

We’ll use a shapefile to crop the raster

To do this we will first get the shapefile of bangalore urban district

ind_shp <- rgdal::readOGR("~/Downloads/IND_adm/IND_adm3.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/actify/Downloads/IND_adm/IND_adm3.shp", layer: "IND_adm3"
## with 2299 features
## It has 13 fields
## Integer64 fields read as strings:  ID_0 ID_1 ID_2 ID_3
DT::datatable(ind_shp@data)
blr_shp <- ind_shp[tolower(ind_shp$NAME_3) == "bangalore",]
plot(blr_shp, 
     border = "blue", 
     main = "Map of Bangalore Urban")

Cropping with shapefile

blr_shp <- spTransform(blr_shp, CRS(proj4string(rasterBrick)))

rasterBrick_cropped <- crop(rasterBrick, blr_shp)

plot(rasterBrick_cropped)

Aggregating and Disaggregating

  • Aggregation is the creation of a lower-resolution raster by grouping rectangular sets of cells in the original raster into individual, larger cells in the new raster.
raster1 <- rasterBrick_cropped[[2]][1:100,1:100, drop = FALSE]

raster1_agg2 <- aggregate(raster1, fact = 2)
raster1_agg4 <- aggregate(raster1, fact = 4)
raster1_agg8 <- aggregate(raster1, fact =8)
raster1_agg16 <- aggregate(raster1, fact = 16)
# raster1 <- raster1[1:100,1:100]

plot(raster1, main = "Original Image")

plot(raster1_agg2, main = "2x2 aggregated")

plot(raster1_agg4, main = "4x4 aggregated" )

plot(raster1_agg8, main = "8x8 aggregated")

plot(raster1_agg16, main = "16x16 aggregated" )


To explain this more clearly, let’s make a few rasters ourselves

set.seed(100)
r <- raster(ncol=10, nrow=10)
r[] <- runif(ncell(r)) * 10
plot(r)
text(r, digits =2 )

r_agg <- aggregate(r, fact = 2)
plot(r_agg)
text(r_agg, digits= 2)

r_disagg <- disaggregate(r_agg, fact = 2)

plot(r_disagg)
text(r_disagg, digits = 2)

Raster Resampling

  • Raster resampling is used to transfer the values of a given raster from its own grid to a different grid.

  • This is done using the resample function

We will first read rasters of different resolutions

raster_B8 <- raster(x = "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1/LC08_L1TP_144051_20181020_20181031_01_T1_B8.TIF")
raster_B8
## class       : RasterLayer 
## dimensions  : 15421, 15101, 232872521  (nrow, ncol, ncell)
## resolution  : 15, 15  (x, y)
## extent      : 632992.5, 859507.5, 1323592, 1554908  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1/LC08_L1TP_144051_20181020_20181031_01_T1_B8.TIF 
## names       : LC08_L1TP_144051_20181020_20181031_01_T1_B8 
## values      : 0, 65535  (min, max)
raster_B2 <- raster(x = "/Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1/LC08_L1TP_144051_20181020_20181031_01_T1_B2.TIF")
raster_B2
## class       : RasterLayer 
## dimensions  : 7711, 7551, 58225761  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 632985, 859515, 1323585, 1554915  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : /Users/actify/Documents/ABG_Training/LC08_L1TP_144051_20181020_20181031_01_T1/LC08_L1TP_144051_20181020_20181031_01_T1_B2.TIF 
## names       : LC08_L1TP_144051_20181020_20181031_01_T1_B2 
## values      : 0, 65535  (min, max)

Next, we will resample raster_B8 to the same rsolution as the raster_B2. There are two methods with which we can do resampling:

  1. “ngb”: Nearest-neighbor; assigns the value of the nearest cell
  2. “bilinear”: Bilinear interpolation; assigns a weighted average of the four nearest cells (the default).

Note : The nearest-neighbor interpolation retains the original values, only transferring them to the new grid according to proximity, while resampling methods that use weighted averages from several cells (such as bilinear interpolation) inevitably involve modification of the original values thus producing a smoother image.

Therefore, in cases where the original values should be preserved, the nearest-neighbor resampling method is recommended.

resampled_raster_B8 <- resample(x = raster_B8, raster_B2, mthod = "ngb")
resampled_raster_B8
## class       : RasterLayer 
## dimensions  : 7711, 7551, 58225761  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 632985, 859515, 1323585, 1554915  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=43 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : LC08_L1TP_144051_20181020_20181031_01_T1_B8 
## values      : 0, 65535  (min, max)